UNCLASSIFIED:  Distribution  Statement  A.  Approved  for  public  release.  #24515 


CONTINUUM  MECHANICS  BASED  BI  LINEAR  SHEAR  DEEORMABLE 
SHELL  ELEMENT  USING  ABSOLUTE  NODAL  COORDINATE 

EORMULATION 


Hiroki  Yamashita 

Department  of  Mechanical  and  Industrial  Engineering 
The  University  of  Iowa 
2312  Seamans  Center 
Iowa  City,  lA  52242 


Antti  I.  Valkeapdd 

Department  of  Mechanical  Engineering 
Lappeenranta  University  of  Technology 
Skinnarilankatu  34, 

53850  Lappeenranta,  Einland 


Paramsothy  Jayakumar 
US  Army  RDECOM  TARDEC 
6501  E.  II  Mile  Road 
Warren,  MI  48397-5000 


Hiroyuki  Sugiyama 

Department  of  Mechanical  and  Industrial  Engineering 
The  University  of  Iowa 
24I6C  Seamans  Center 
Iowa  City,  lA  52242 
hiroyuki- sugiyama  @  uiowa.  edu 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  0MB  control  number. 

1 .  REPORT  DATE  2.  REPORT  TYPE 

07  MAR  2014  Technical  Report 

4.  TITLE  AND  SUBTITLE 

CONTINUUM  MECHANICS  BASED  BI  LINEAR  SHEAR 
DEFORMABLE  SHELL  ELEMENT  USING  ABSOLUTE  Nt 
COORDINATE  FORMULATION 

6.  AUTHOR(S) 

Hiroki  Yamashita;  Antti  Valkeapaa;  Paramsothy  Jayakumar: 

Sugiyama 

7.  PEREORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Department  of  Mechanical  and  Industrial  Engineering, The  Ui 
of  Iowa, 2312  Seamans  Center,Iowa  City,IA,52242 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  TARDEC,  6501  East  Eleven  Mile  Rd,  Warren,  Mi, 

48397-5000 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

For  ASME  Journal  of  Computational  Nonlinear  Dynamics 

14.  ABSTRACT 

In  this  investigation,  a  continuum  mechanics  based  bi-linear  shear  deformable  shell  element  is  developed 
using  the  absolute  nodal  coordinate  formulation  for  the  large  deformation  analysis  of  multibody  shell 
structures.  The  element  consists  of  four  nodes,  each  of  which  has  the  global  position  coordinates  and  the 
gradient  coordinates  along  the  thickness  introduced  for  describing  the  orientation  and  deformation  of  the 
cross  section  of  the  shell  element.  The  global  position  field  on  the  mid-surface  and  the  position  vector 
gradient  at  a  material  point  in  the  element  are  interpolated  by  bi-linear  polynomials.  The  continuum 
mechanics  approach  is  used  to  formulate  the  generalized  elastic  forces,  allowing  for  the  consideration  of 
nonlinear  constitutive  models  in  a  straightforward  manner.  The  element  locking  exhibited  in  this  type  of 
element  can  be  eliminated  using  the  assumed  natural  strain  (ANS)  and  enhanced  assumed  strain  (EAS) 
approaches.  In  particular,  the  combined  ANS  and  EAS  approach  is  introduced  to  alleviate  the  thickness 
locking  arising  from  the  erroneous  transverse  normal  strain  distribution.  Several  numerical  examples  are 
presented  in  order  to  demonstrate  the  accuracy  and  the  rate  of  convergence  of  numerical  solutions 
obtained  by  the  continuum  mechanics  based  bi-linear  shear  deformable  ANCF  shell  element  proposed  in 
this  investigation. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF 

OF  ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Public 

Release 

42 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


UNCLASSIFIED:  Distribution  Statement  A.  Approved  for  public  release.  #24515 


ABSTRACT 


In  this  investigation,  a  continuum  mechanics  based  bi-linear  shear  deformable  shell  element  is 
developed  using  the  absolute  nodal  coordinate  formulation  for  the  large  deformation  analysis  of 
multibody  shell  structures.  The  element  consists  of  four  nodes,  each  of  which  has  the  global 
position  coordinates  and  the  gradient  coordinates  along  the  thickness  introduced  for  describing 
the  orientation  and  deformation  of  the  cross  section  of  the  shell  element.  The  global  position 
field  on  the  mid-surface  and  the  position  vector  gradient  at  a  material  point  in  the  element  are 
interpolated  by  bi-linear  polynomials.  The  continuum  mechanics  approach  is  used  to  formulate 
the  generalized  elastic  forces,  allowing  for  the  consideration  of  nonlinear  constitutive  models  in  a 
straightforward  manner.  The  element  locking  exhibited  in  this  type  of  element  can  be  eliminated 
using  the  assumed  natural  strain  (ANS)  and  enhanced  assumed  strain  (EAS)  approaches.  In 
particular,  the  combined  ANS  and  EAS  approach  is  introduced  to  alleviate  the  thickness  locking 
arising  from  the  erroneous  transverse  normal  strain  distribution.  Several  numerical  examples  are 
presented  in  order  to  demonstrate  the  accuracy  and  the  rate  of  convergence  of  numerical 
solutions  obtained  by  the  continuum  mechanics  based  bi-linear  shear  deformable  ANCF  shell 
element  proposed  in  this  investigation. 
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1.  INTRODUCTION 

The  rectangular  plate/shell  elements  of  the  absolute  nodal  coordinate  formulation  (ANCF)  [1] 
can  be  classified  into  the  fully  parameterized  shear  deformable  element  [2,3]  and  the  gradient 
deficient  thin  plate  element  [4,5].  While  the  fully  parameterized  element  leads  to  a  general 
motion  description  that  accounts  for  complex  coupled  deformation  modes  of  the  plate/shell 
element,  use  of  higher  order  polynomials  and  the  coupled  deformation  modes  exhibited  in  this 
type  of  element  causes  severe  element  locking  and  special  care  needs  to  be  exercised  to  alleviate 
the  locking  [6].  The  gradient  deficient  thin  plate  element,  on  the  other  hand,  is  developed  by 
removing  the  position  vector  gradients  along  the  thickness  {dr ! dz),  and  the  global  displacement 
field  on  the  mid-surface  can  be  uniquely  parameterized  by  the  global  position  vector  and  the  two 
gradient  vectors  dr  !  dx  and  dr  !  dy  which  are  both  tangent  to  the  surface.  By  doing  so,  the  cross 
section  is  assumed  to  be  rigid  and  the  elastic  forces  are  derived  using  an  in-plane  stress 
assumption  with  Kirchhoff-Love  plate  theory.  This  leads  to  the  non-conforming  plate  element  in 
which  the  inter-element  continuity  is  guaranteed,  while  the  conforming  thin  plate  element  can  be 
obtained  by  introducing  the  additional  nodal  coordinates  d\ !  dxdy  with  bi-cubic  Hermite 
polynomials  [4].  Despite  the  fact  that  the  gradient  deficient  thin  ANCF  plate  elements  have 
proven  to  be  successful  in  solving  challenging  engineering  problems  that  involve  large 
deformable  thin  plate  and  shell  structures,  consideration  of  general  nonlinear  constitutive  models 
requires  special  formulations  and  implementation  due  to  the  in-plane  stress  assumption. 

In  recent  years,  an  ANCF  parameterization  which  does  not  include  the  position  vector 
gradients  tangent  to  the  beam  centerline  and  the  mid-surface  are  investigated  for  shear 
deformable  ANCF  beam  and  plate  elements.  For  beam  problems,  it  is  shown  in  the  literature  [7] 
that  the  elimination  of  the  tangential  slope  vector  along  the  beam  centerline  leads  to  accurate 
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elastic  force  descriptions  due  to  weaker  polynomial  coupling  between  the  displacement  and 
gradient  fields,  and  use  of  the  general  continuum  mechanics  approach  for  the  elastic  forces 
calculation  allows  for  modeling  curved  beam  structures  in  a  straightforward  manner.  The 
assumed  displacement  field  on  the  beam  centerline  does  not  involve  slope  coordinates,  while  the 
transverse  slopes  are  employed  to  parameterize  the  orientation  and  deformation  of  the  beam 
cross  section.  For  plate  elements,  on  the  other  hand,  the  position  vector  gradient  along  the 
thickness  can  be  used  to  describe  the  orientation  and  deformation  of  the  cross  section  of  the  plate 
element,  and  this  parameterization  leads  to  a  shear  deformable  element  that  can  capture  the 
thickness  deformation  using  the  transverse  slope  vectors.  That  is,  use  of  such  an  element 
parameterization  leads  to  shear  deformable  plate  elements  while  reducing  the  number  of 
coordinates  per  node  as  discussed  in  the  literature  [8],  in  which  a  selective  reduced  integration  is 
used  to  alleviate  the  transverse  shear  locking  of  the  bi-linear  shear  deformable  ANCF  plate 
element.  The  element  locking  problem  exhibited  in  the  bi-linear  shear  deformable  ANCF  plate 
element  is  further  investigated  in  the  literature  [9],  in  which  the  assumed  natural  strain  (ANS) 
and  the  enhanced  assumed  strain  (EAS)  approaches  are  suggested  to  alleviate  the  transverse 
shear  locking,  thickness  locking,  and  the  in-plane  normal/shear  locking. 

The  elastic  forces  of  the  bi-linear  shear  deformable  ANCF  plate  element  in  the  literatures 
[8,9]  are  formulated  based  on  the  elastic  mid-plane  approach,  in  which  the  six  strain  components 
are  evaluated  on  the  mid-plane  using  Green-Lagrange  strains  and  the  approximated  curvature 
expression  is  used  to  define  the  plate  bending  and  twisting  deformation.  In  other  words,  the 
strain  distribution  along  the  thickness  is  assumed  to  be  constant  in  this  approach,  allowing  for  the 
evaluation  of  the  elastic  forces  as  an  area  element.  Use  of  this  element  is  limited  to  moderately 
thick  flat  plate  problems  due  to  the  simplifying  assumptions  made  in  this  model.  Furthermore, 
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consideration  of  general  nonlinear  constitutive  models  requires  special  formulations  and 
implementation  as  in  the  thin  ANCF  plate  elements.  Most  challenging  engineering  problems  that 
involves  large  deformable  shell  structures,  on  the  other  hand,  call  for  accurate  modeling  of  shell 
structures  with  various  material  nonlinearities,  some  of  which  include  the  modeling  of  rubber 
tires  for  vehicle/terrain  interaction;  wind  turbine  rotor  blades  made  of  composite  materials; 
biomechanics  applications  with  biomaterials;  and  many  others.  To  achieve  this  goal,  elastic 
forces  of  the  shell  element  need  to  be  developed  as  a  continuum  volume  by  considering  the  stress 
variation  along  the  thickness  using  a  general  continuum  mechanics  approach,  while  the  thickness 
locking  can  have  an  adverse  impact  on  the  accuracy  as  has  been  address  for  solid  shell  elements 
[10-13].  It  is,  therefore,  the  objective  of  this  investigation  to  develop  the  continuum  mechanics 
based  bi-linear  shear  deformable  shell  element  using  the  absolute  nodal  coordinate  formulation. 

The  paper  that  discusses  the  development  of  the  bi-linear  shear  deformable  ANCF  shell 
element  using  the  continuum  mechanics  approach  is  organized  as  follows.  In  Section  2,  the 
kinematics  and  parameterization  of  the  bi-linear  ANCF  shell  element  are  presented,  while  the 
elastic  forces  based  on  the  general  continuum  mechanics  approach  and  remedies  for  the  element 
locking  exhibited  in  the  element  are  discussed  in  Section  3.  In  Section  4,  the  equations  of  motion 
and  the  solution  procedure  are  discussed.  Numerical  results  and  comparison  with  existing 
plate/shell  elements  are  presented  in  order  to  demonstrate  the  effect  of  the  element  locking  on  the 
accuracy  for  the  continuum  mechanics  based  bi-linear  shear  deformable  ANCF  shell  element  in 
Section  5.  Summary  and  conclusions  drawn  from  this  study  are  presented  in  Section  6. 

2.  KINEMATICS  OF  BI-LINEAR  SHEAR  DEFORMABLE  ANCF  SHELL  ELEMENT 

As  shown  in  Fig.  1,  the  global  position  vector  r'  of  a  material  point  x'  =[v‘  y'  in  shell 
element  i  is  defined  as 
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OZ 


(1) 


where  ,y)  is  the  global  position  vector  on  the  mid-surface  and  drl^{x' ,y')l dz'  is  the 

transverse  gradient  vector  used  to  describe  the  orientation  and  deformation  of  the  infinitesimal 
volume  in  the  element.  The  preceding  global  displacement  field  is  interpolated  using  the  bi-linear 
polynomials  as  follows: 


r‘  =aQ  -hfljx'  -t-fljy  +a.j,x'y  +z‘{a^  +a^x'  +a^y  +a-,x'y) 


(2) 


from  which,  one  can  interpolate  both  position  vector  on  the  mid-surface  and  the  transverse 
gradient  vector  using  the  same  bi-linear  shape  function  matrix  as  follows: 


<  (X' ,  v' ) = s;.  (y ,  y  )e;. ,  ^  ,  y  > = si  (y ,  y 


(3) 


where  S;„=[5;i  5'I  5‘I  5'l]  and 


where  =x'  !  (i‘  /  2)  and  rj'  =  y'  !  (w‘  /  2) .  t  and  w'  are  lengths  along  the  element  x‘  and  y‘ 
axes,  respectively.  In  Eq.  3,  the  vectors  and  e'^  represent  the  element  nodal  coordinates 
associated  with  the  global  position  vector  on  the  mid-surface  and  the  transverse  gradient  vector. 
That  is,  for  node  k  of  element  i,  one  has  and  dz’ .  It  is  important  to  notice  here 

that  the  assumed  displacement  field  defined  on  the  mid-surface  in  Eq.  3  does  not  involve  any 

gradient  coordinates,  while  the  orientation  and  deformation  of  the  infinitesimal  volume  at  the 
material  point  in  the  shell  element  is  parameterized  by  the  transverse  gradient  coordinates  only. 
Substitution  of  Eq.  3  into  Eq.  1  leads  to  the  following  general  expression  for  the  global  position 
vector  used  for  the  absolute  nodal  coordinate  formulation: 
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r'(x'',y,z')  =  S'(x'',y,z')e''  (5) 

where  the  shape  function  matrix  S'  and  the  element  nodal  coordinate  vector  e'  are,  respectively, 
defined  as 


S'=[SL  j'SL],  e'=[(e:,f  (e;fr  (6) 

3.  CONTINUUM  MECHANICS  BASED  SHEAR  DEEORMABLE  ANCF  SHELL 
ELEMENT 


3.1  Generalized  Elastic  Forces  using  Continuum  Mechanics  Approach 


In  the  continuum  mechanics  approach,  the  elastic  forces  of  the  shell  element  are  evaluated  as  a 
continuum  volume  and  the  Green-Lagrange  strain  tensor  E  at  an  arbitrary  material  point  in 
element  i  is  defined  as  follows: 


E''  =-((F'yF'-l) 


(7) 


where  F'  is  the  displacement  gradient  tensor  defined  by 


F'  = 


dr'  _dr‘ 
dX‘  ~  dx‘ 


i  f  csvi  ^ 


(8) 


In  the  preceding  equation,  J'  =dr‘/dx‘  and  J'  =5X'/5x' ,  where  the  vector  X'  represents  the 
global  position  vector  of  element  i  at  an  arbitrary  reference  configuration.  Substitution  of  Eq.  8 
into  Eq.  7  leads  to 

E'=(J')-^E'(J')-'  (9) 

where  E'  is  the  covariant  strain  tensor  defined  by 


E'=i((j'yj'-(j'yj') 


(10) 
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The  transformation  (push-forward  operation)  of  the  covariant  strain  tensor  E'  given  by  Eq.  9  can 
be  re-expressed  in  a  vector  form  by  introducing  the  engineering  covariant  strain  vector  e'  as 

2'=(T)-^£'  (11) 

where  the  vector  e  in  Eq.  11  is  the  engineering  strain  vector  at  the  deformed  configuration 
defined  as 


e' =[s‘  s'  /  s'  /  /f 

>-  XX  yy  /  4!^/  zz  /  xz  /  vz  -■ 


(12) 


and  the  engineering  covariant  strain  vector  is  defined  as 

£''=[4  4  4  4  4  fyzf 

The  transformation  matrix  T'  in  Eq.  1 1  can  be  expressed  explicitly  as 


(13) 


rp)  _ 


iJnf 

iJnf 

2J 

i  ji 
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(14) 


and  7'^  is  the  element  in  the  a-th  column  and  b-i\\  row  of  matrix  J'  which  is  constant  in  time. 


The  generalized  elastic  forces  can  then  be  obtained  using  the  virtual  work  as  follows: 


« =L 


c'dV^ 


(15) 


where  a'  is  a  vector  of  the  second  Piola-Kirchhoff  stresses  and  dV^  is  the  infinitesimal  volume 

at  the  reference  configuration  of  element  i.  It  is  important  to  notice  here  that  the  element  elastic 
forces  are  evaluated  as  a  continuum  volume,  and  the  stress  vector  for  the  shell  element  can  be 
obtained  with  various  nonlinear  constitutive  models  for  large  deformation  problems  without  ad 
hoc  procedures. 
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3.2  Element  Locking  for  Transverse  Shear  and  In-Plane  Shear/Normal  Strains 

As  it  has  been  addressed  in  literatures  of  shell  element  formulations  [10-20],  the  bi-linear 
quadrilateral  shell  element  suffer  from  the  transverse  shear  locking  and  the  in-plane  shear/normal 
locking.  The  transverse  shear  locking  can  be  eliminated  using  the  assumed  natural  strain  (ANS) 
approach  proposed  by  Bathe  and  Dvorkin  [16,17].  In  this  approach,  the  covariant  shear  strain 
components  are  interpolated  using  those  evaluated  at  the  sampling  points  A,  B,C  and  D  shown  in 
Fig.  2  as  follows: 

(16) 

where  f  ® and  f  ®  are  compatible  covariant  transverse  shear  strains  at  the  sampling 

points. 

The  parasitic  in-plane  shear  under  pure  bending  loads  is  a  typical  locking  problem 
exhibited  in  the  bi-linear  quadrilateral  element  [21],  and  the  compatible  in-plane  strains 

and  y^  )  obtained  by  the  assumed  displacement  field  can  be  enhanced  by  introducing  the 

enhanced  assumed  strains  (EAS)  as  [14,15] 

£  =  £‘^+£“^  (17) 

where  e"  indicates  the  compatible  strain  vector  and  the  strain  vector  £^''^  is  defined  by 

£“"(^)  =  G(§)a  (18) 

In  the  preceding  equation,  a  is  a  vector  of  internal  parameters  introduced  to  define  the  enhanced 
in-plane  strain  field  and  the  matrix  G(^)  can  be  defined  as  [14,15] 
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=  (19) 

where  J(^)  and  Jg  are  displacement  gradient  matrices  at  the  reference  configuration  evaluated 
at  the  Gaussian  integration  point  ^  and  at  the  center  of  element  ( ^  =  0 ),  respectively.  ^  is  a 
vector  of  the  element  coordinates  in  the  parametric  domain  and  Tg  is  the  constant  transformation 
matrix  evaluated  at  the  center  of  element  as  shown  in  Eq.  14.  The  matrix  M(^)  defines 
polynomials  for  the  enhancement  of  the  in-plane  strain  field  in  the  parametric  domain.  For 
example,  the  linear  distribution  of  in-plane  strains  (e^,  and  )  requires  introducing  the 

following  interpolation  matrix  M(^) : 

^000 
0^00 
0  0  ^  ri 

0  0  0  0 
0  0  0  0 
0  0  0  0 

where  the  additional  four  internal  FAS  parameters  associated  with  the  four  columns  of  the  matrix 
M(^)  are  introduced  in  this  model.  Using  Eq.  19,  the  enhanced  covariant  strains  are  then  pushed 
forward  to  those  at  the  deformed  configuration  in  the  physical  domain.  It  is  important  to  notice 
here  that  the  matrix  M(^)  needs  to  satisfy  the  following  condition  [14]: 

jM(§)rf^  =  0  (21) 

such  that  the  orthogonality  condition  between  the  assumed  stress  and  strain  is  satisfied  as 

f  o-£“^^/Ug=0  (22) 

JVo 
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Using  the  preceding  condition,  the  unknown  assumed  stress  term  that  appears  in  Hu-Washizu 
mixed  variational  principle  vanishes  and  one  can  obtain  the  generalized  elastic  force  vector  as 
follows  [14]: 


ds‘ 


dV‘ 


(23) 


where  W  is  an  elastic  energy  function.  For  example,  the  elastic  energy  function  W  of  the 
compressible  neo-Hookean  material  model  is  defined  as  follows  [22]: 


W  =  ^(tr(C)-3)-//ln/  +  ^(lny)" 


A 


(24) 


where  //  and  A  are  the  Lame  constants;  C  =  F^F  is  the  right  Cauchy-Green  deformation  tensor; 
and  J  =  det(F)  =  ^Jdet(C) .  It  is  important  to  notice  here  that  the  right  Cauchy-Green  deformation 
tensor  that  accounts  for  the  enhanced  assumed  strain  modification  is  defined  as  follows: 


C  =  2(E"+E'^^)+I 


(25) 


3.3  Thickness  Locking 

Use  of  the  transverse  gradients  in  the  ANCF  shell  element  introduces  thickness  deformation 
modes  and  the  locking  associated  with  the  transverse  normal  strain  is  exhibited  in  the  bi-linear 
shear  deformable  ANCF  shell  element.  The  solid  shell  elements  which  consists  of  layers  of 
translational  nodal  coordinates  at  the  top  and  bottom  surfaces  of  the  element  [10-13]  and  shell 
elements  parameterized  by  extensible  directors  [18-20]  suffer  from  the  similar  thickness  locking 
resulting  from  the  erroneous  transverse  normal  strain  distribution.  The  use  of  the  assumed  natural 
strain  (ANS)  approach  is  proposed  in  the  literature  [18]  to  alleviate  the  thickness  locking  of  the 
shell  element  modeled  by  extensible  directors  and  the  transverse  normal  strain  at  a  material  point 
in  the  element  is  approximated  as  follows: 
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_  ^ANS  1  ^ANS  2  ^ANS  3  ^ANS  4 
“  "^1  ^zz  ^3  ^zz  "^4  ^zz 


(26) 


where  indicates  the  compatible  transverse  normal  strain  at  node  k  as  shown  in  Fig.  2  and 


is  the  shape  function  associated  with  it.  This  approach  is  applied  to  the  bi-linear  shear 

deformable  ANCF  flat  plate  element  formulated  with  the  elastic  mid-plane  approach  [9]  to 
alleviate  the  thickness  locking  problem.  Since  the  strain  distribution  along  the  thickens  is 
assumed  to  be  constant  and  the  stress  is  evaluated  on  the  mid-surface,  use  of  the  assumed  natural 
strain  approach  can  alleviate  the  thickness  locking  effectively  if  the  compatible  transverse 
normal  strains  at  four  corners  on  the  mid-surface  are  accurate  enough.  However,  in  the  case  of 
the  continuum  mechanics  approach,  the  elastic  forces  are  evaluated  as  a  continuum  volume  and 
the  transverse  normal  strain  are  evaluated  at  the  integration  points  which  do  not  lie  on  the  mid¬ 
surface,  thus  better  approximation  of  the  transverse  normal  strain  is  needed  to  alleviate  the 
thickness  locking  for  the  continuum  mechanics  based  ANCF  shell  element. 

It  has  been  shown  that  the  use  of  the  enhanced  assumed  strain  (EAS)  approach  leads  to 
the  better  approximation  of  the  transverse  normal  strain  distribution  and  the  thickness  locking 
can  be  alleviated  [11,19].  To  this  end,  the  additional  internal  EAS  parameters  associated  with  the 
enhanced  transverse  normal  strain  are  added,  and  the  interpolation  matrix  M(^)  is  augmented 
associated  with  it.  For  example,  if  the  linear  transverse  normal  strain  term  is  introduced  as  a 
function  of  the  thickness  coordinate  z,  the  interpolation  matrix  M(^)  given  in  Eq.  19  is  modified 
as 
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^  0  0  0  0 

0  7  0  0  0 

0  0  ^  7  0 

0  0  0  0  C 

0  0  0  0  0 

0  0  0  0  0 

where  =  z/(h/ 2)  .  It  is  important  to  notice  here  that  there  are  no  strain  enhancement 
associated  with  the  transverse  shear  strains  and  defined  in  the  fifth  and  sixth  rows.  While 
there  are  zero  sub-matrix  associated  with  these  strains  in  the  matrix,  the  matrix  G(^)  define  by 
Eq.  19  becomes  a  full  matrix  since  the  matrix  Tq“^  is  defined  as  a  6x6  full  matrix  evaluated  at 
an  arbitrary  reference  configuration.  In  the  literature  [19],  the  transformation  matrix  G(^) 


associated  with  the  enhanced  in-plane  and  thickness  strain  terms  are  assumed  to  be  decoupled 
and  the  approximated  transformation  matrix  is  used.  In  this  investigation,  the  exact 
transformation  matrix  defined  for  the  general  three-dimensional  stress  state  is  considered. 

To  further  improve  the  transverse  normal  strain  distribution  in  the  continuum  mechanics 
based  element,  the  compatible  strain  for  the  transverse  normal  strain  is  modified  by  the  assumed 
natural  strain  approach  defined  by  Eq.  26  as  follows: 


(28) 


The  similar  approach  is  employed  for  solid  shell  elements  in  literature  [11].  In  this  investigation, 
the  combined  ANS  and  EAS  approach  for  the  thickness  locking  is  applied  to  the  bi-linear  shear 
deformable  ANCF  shell  element.  With  this  approach,  the  transverse  normal  strain  improved  by 
the  ANS  approach  is  further  enhanced  by  the  enhanced  assumed  strain  approach.  Eor  application 
of  the  approach  to  the  continuum  mechanics  based  ANCE  shell  element,  the  covariant  strain 


components  of  the  transverse  shear  and  transverse  normal  strains  are  interpolated  in  the  natural 
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coordinate  domain  first  and  then  the  covariant  strain  vector  given  in  Eq.  13  is  replaced  with  the 
following  strain  vector: 


The  Green-Lagrange  strains  associated  with  the  preceding  covariant  strains  are  evaluated  using 
the  transformation  defined  by  Eq.  11.  Having  obtained  the  modified  compatible  strains,  the 
enhanced  assumed  strains  defined  by  Eq.  18  that  considers  the  strain  enhancement  of  the  in¬ 
plane  shear/normal  as  well  as  the  transverse  normal  strains  are  added  to  them  with  the 
interpolation  matrix  M(§)  given  by  Eq.  27.  This  leads  to  a  systematic  derivation  of  the 


generalized  elastic  forces  of  the  locking-free  continuum  mechanics  based  shear  deformable 
ANCF  shell  element. 


4.  EQUATIONS  OF  MOTION  AND  SOLUTION  PROCEDURES 

Using  the  principle  of  virtual  work  in  dynamics,  the  equations  of  motion  of  the  shear  deformable 
ANCF  shell  element  i  can  be  expressed  as 

M'e'=QUe'',a')+Q:(e'',e'',0  (30) 

where  the  vectors  and  are,  respectively,  the  element  elastic  and  external  force  vectors; 
and  the  matrix  M'  is  the  constant  element  mass  matrix  defined  by 

M'  =  f  ^(S'/SW'  (31) 

where  is  the  material  density  at  the  reference  configuration.  It  is  important  to  notice  here  that 
the  generalized  elastic  forces  are  defined  as  a  function  of  the  nodal  coordinates  and  the  internal 
FAS  parameters  a' .  The  parameter  a'  are  determined  by  solving  the  following  equation  for  each 
element: 
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h'(e',a')=f 

JV( 


dV‘  =  0, 


i  =  l,---,ne 


The  nodal  coordinates  e'  and  the  EAS  internal  parameters  a'  are  determined  at  every  time  step 
such  that  both  Eqs.  30  and  32  are  satisfied.  Using  the  time  integration  scheme  in  Newmark-|3 
family,  the  nodal  acceleration  and  velocity  vectors  at  time  step  are  discretized  in  time  as 


®n+l 


1  ,  ,1.1  .... 
3,  2  nu  ®n  L  o 


Ph  "  IP 


e„,  =  ^  (e„,  -  e. )  -  (^  - 1)  e.  -  -  l)e„ 

Using  the  preceding  equations,  the  equations  of  motion  given  by  Eq.  30  can  be  expressed  as  a 
function  of  unknown  variables  and  a.\^^  as  follows: 

f‘'(eLp«L)  =  0  (34) 

In  other  words,  one  needs  to  seek  solutions  to  the  following  equations  at  time  step  : 


f '  te'  a!  1 

sTe”'  a'  )= 


Using  Newton-Raphson  method,  Newton  difference  for  and  can  be  obtained  by  solving 
the  following  matrix  equation: 

[dl'/deL.,  SfVSa'  ,irAei.,l  frl 


Sf/SeL,  dl‘/dal 
ah'/aei.,  ah'/aa;. 


Using  the  preceding  equation,  can  be  obtained  as  a  solution  to  the  following  equation: 


■'  5h' 

5f' 

Ae'  ,  =  f'  , 

Sell 

8e' 

^^n+l  J 

n+l  n+l  ^  i 

The  preceding  equation  is  iteratively  solved  until  g'(e^^j,a^^j)  =  0  defined  by  Eq.  35  is  satisfied. 
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The  internal  EAS  parameters  are  determined  in  each  iteration  step  using  the  following  update 
scheme: 


a;,^=a;-(H‘rh‘(e‘,a;) 

where  Jacobian  matrix  H'  in  the  preceding  equation  is  defined  by: 

ah'  f  aV'(£"+s“^) 


(38) 


H‘  =  — =f  (G'Z 
da‘ 


-G‘dV: 


(39) 


where  the  matrix  G'  is  defined  by  Eq.  19.  The  term  d^W  /  leads  to  a  material  moduli  matrix 
obtained  by  the  second  derivative  of  the  elastic  energy  function.  That  is,  for  a  linear  elasticity 
material  model,  this  matrix  is  constant.  However,  for  general  nonlinear  material  models,  the 
material  moduli  matrix  is  no  longer  constant  and  the  iterative  solution  procedure  is  required  to 
determine  the  internal  EAS  parameters  a'  for  each  element. 

5.  NUMERICAL  RESULTS 


In  this  section,  several  numerical  examples  are  presented  in  order  to  demonstrate  the 
performance  of  the  continuum  mechanics  based  bi-linear  shear  deformable  ANCE  shell  element 
developed  in  this  investigation.  The  effect  of  the  assumed  natural  strain  (ANS)  and  enhanced 
assumed  strain  (EAS)  approaches  on  the  element  accuracy  is  also  discussed. 

5.1  Cantilevered  Plate  and  Shell  Subjected  to  a  Point  Force 

In  the  first  problem,  a  rectangular  cantilevered  plate  subjected  to  a  vertical  point  force  at  one  of 
the  corners  of  the  plate  is  considered  as  shown  in  Fig.  3.  The  length,  width  and  thickness  of  the 
plate  are  assumed  to  be  1.0  m,  1.0  m  and  0.01  m.  The  Young’s  modulus  and  Poisson’s  ratio  are 
assumed  to  be  2.1x10*  Pa  and  0.3,  respectively.  As  shown  in  Fig.  3,  the  plate  is  subjected  to 
large  deformation  at  the  static  equilibrium  state.  The  six  models  with  different  strain 
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modifications  discussed  in  Section  3  are  considered  to  demonstrate  the  effect  of  the  EAS  and 
ANS  approaches  on  the  element  convergence  and  accuracy.  These  models  are  summarized  in 
Table  1.  In  Model- 1,  no  strain  modifications  are  made,  while  Model-2  is  used  to  demonstrate  the 
effect  of  the  ANS  approaeh  to  alleviate  the  transverse  shear  locking  for  the  ANCF  shell  element. 
A  eomparison  between  Model-3  and  4  demonstrates  the  effect  of  the  in-plane  shear/normal 
loeking  alleviated  by  the  EAS  approaeh,  while  a  eomparison  between  Model-4,  5  and  6 
demonstrates  the  use  of  the  three  different  approaehes  for  elimination  of  the  thiekness  loeking, 
i.e.,  the  ANS  (Model-4),  the  EAS  (Model-5),  and  the  combined  ANS  and  EAS  (Model-6) 
approaches  as  discussed  in  Section  3.3. 

To  measure  the  accuracy,  error  in  solution  is  defined  by  a  deviation  of  the  vertical 
deflection  at  the  force  application  point  on  the  plate  from  that  of  the  reference  solution  obtained 
by  ANSYS  SHELL181  with  100x100  elements.  The  numerieal  convergenee  of  the  solution  is 
presented  in  Fig.  4  for  each  model.  The  numerieal  result  obtained  by  the  bi-linear  shear 
deformable  ANCF  flat  plate  element  based  on  the  elastie  mid-plane  approach  [9]  is  also 
presented  in  this  figure  for  comparison.  It  is  observed  from  the  comparison  between  Model- 1  and 
Model-2  that  the  accuracy  is  improved  by  alleviating  the  transverse  shear  locking  with  the  ANS 
approach.  However,  the  solution  is  converged  to  the  erroneous  solution.  While  the  use  of  the 
EAS  for  the  in-plane  shear/normal  locking  and  the  ANS  for  the  thickness  locking  slightly 
improve  the  rate  of  convergence,  the  convergent  solution  is  still  incorrect.  It  is  important  to 
notice  that  use  of  Model-6  with  the  combined  ANS  and  EAS  approach  leads  to  the  same  result 
with  that  of  the  loeking-free  ANCF  plate  element  based  on  the  elastic  mid-plane  approaeh.  This 
plate  element  employs  ANS  for  alleviating  the  thickness  loeking.  In  other  words,  the  strain 
modification  used  in  the  elastie  mid-plane  ANCF  shell  element  corresponds  to  Model-4  in  Table 
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1.  However,  Model-4  of  the  continuum  mechanics  based  ANCF  shell  element  leads  to  an 
erroneous  solution  due  to  the  thickness  locking.  While  the  application  of  the  EAS  approach  to 
the  thickness  locking  improves  the  element  performance  as  shown  in  the  result  of  Model-5,  the 
accuracy  is  not  at  satisfactory  level  especially  when  the  small  number  of  elements  is  used.  This 
result  clearly  indicates  severity  of  the  thickness  locking  exhibited  in  the  continuum  mechanics 
based  shear  deformable  ANCF  shell  element  and  the  locking  for  the  continuum  mechanics  based 
ANCF  shell  element  needs  to  be  eliminated  using  the  combined  EAS  and  ANS  approach  as 
demonstrated  by  the  result  of  Model-6. 

In  the  second  problem,  the  flat  plate  is  replaced  with  the  quarter  cylinder  modeled  by  the 
continuum  mechanics  based  bi-linear  ANCE  shell  elements  as  shown  in  Fig.  5,  and  the  radius  of 
curvature  is  assumed  to  be  1.0  m.  The  material  properties,  tip  load,  width  and  height  are  the  same 
as  the  previous  example.  The  deformed  shape  at  the  static  equilibrium  state  is  shown  in  Fig.  5 
and  the  large  deformation  is  exhibited  in  this  problem.  As  in  the  previous  example,  error  in 
solution  is  defined  by  a  deviation  of  the  vertical  deflection  at  the  force  application  point  from 
that  of  the  reference  solution  obtained  by  ANSYS  SHELL181  with  100x100  elements.  The 
numerical  convergence  in  solution  is  presented  in  Fig.  6  for  Models  1  to  6  and  ANSYS 
SHELL181  for  comparison.  The  overall  trend  of  the  numerical  convergence  of  solutions 
obtained  by  different  models  is  similar  to  that  observed  for  the  flat  plate  problem.  The  result 
clearly  indicate  that  use  of  the  combine  EAS  and  ANS  approach  leads  to  the  ideal  rate  of 
convergence  and  accuracy  for  the  shell  structure. 

5.2  Pinched  Semi-Cylindrical  Shell 

A  semi-cylindrical  shell  as  shown  in  Fig.  7  subjected  to  a  pinching  force  of  800  N  at  the  middle 
of  the  free-hanging  circumferential  periphery  is  discussed.  This  problem  is  used  as  a  benchmark 
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problem  for  various  shell  elements  [13,23,24,25]  and  the  numerical  solution  obtained  by  the 
proposed  continuum  mechanics  based  bi-linear  shear  deformable  shell  element  (Model-6)  is 
compared  with  those  in  the  literatures  [13,24,25].  The  length,  radius  and  thickness  are  assumed 
to  be  3.048  m,  1.016  m  and  0.03  m,  respectively.  Young’s  modulus  and  Poisson’s  ratio  are 
assumed  to  be  2.0685x10^  Pa  and  0.3,  respectively  [13,23,24,25].  A  quarter  cylinder  model  with 
the  symmetric  boundary  condition  is  used  in  the  analysis.  The  deformed  shape  of  the  semi- 
cylindrical  shell  is  shown  in  Fig.  7,  and  the  nonlinear  load-deflection  curve  shown  in  Fig.  8  agree 
well  with  those  presented  in  the  literatures  [13,24,25].  Use  of  the  16x16  elements  leads  to  the 
convergent  solution  in  the  case  of  the  continuum  mechanics  based  bi-linear  shear  deformable 
ANCF  shell  element. 

5.3  Slit  Annular  Plate  Subjected  to  Lifting  Force 

In  the  next  example,  a  slit  annular  plate  subjected  to  a  lifting  line  force  is  considered  as  shown  in 
Fig.  9.  This  problem  is  also  widely  used  as  a  benchmark  problem  for  shell  elements  [23].  The 
inside  and  outside  radii  of  the  annular  plate  are  assumed  to  be  6.0  m  and  10.0  m,  respectively. 
Young’s  modulus  and  Poisson’s  ratio  are  assumed  to  be  2.1x10’  Pa  and  v  =  0.0,  respectively 
[23].  One  end  of  the  slit  is  fully  clamped,  while  the  other  end  of  the  slit  is  subjected  to  the  line 
force  of  50  N.  The  load-deflection  curves  at  Point  A  and  B  shown  in  Fig.  10  are  compared  with 
those  obtained  using  ANSYS  SHELL181  with  10x60  elements  as  a  reference  solution  and  the 
results  are  in  good  agreement  with  those  obtained  by  the  continuum  mechanics  based  bi-linear 
shear  deformable  ANCF  shell  element. 

5.4  Natural  Frequencies  of  Square  Plate 

The  eignefrequency  analysis  of  a  square  plate  with  free  boundaries  studied  in  the  literatures 
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[6,26]  is  performed  in  this  example.  The  dimensionless  natural  frequencies  defined  by 

Q  =  co/a>Q  with  cOq  =  pht  are  used,  where  D,  p,  h  and  I  are  D  =  Eh^  l\2{\-v^),  the 

material  density,  height  and  length  of  the  plate,  respectively.  The  first  ten  dimensionless 
eigenfrequencies  are  summarized  in  Table  2  for  Model-6  that  employs  the  combined  ANS  and 
EAS  approach  for  the  thickness  locking.  Use  of  the  proposed  low-order  continuum  mechanics 
based  shear  deformable  shell  element  leads  to  better  rate  of  convergence  than  the  higher-order 
fully  parameterized  ANCF  element  presented  in  literatures  [6,26]. 

5.5  Quarter  Cylinder  Pendulum 

In  the  last  problem,  the  nonlinear  dynamic  analysis  of  a  shell  structure  with  a  nonlinear  material 
model  is  considered.  A  quarter  cylinder  that  has  the  same  dimension  as  the  one  in  Section  5.1  is 
used,  and  a  compressible  neo-Hookean  material  model  is  used  to  demonstrate  the  use  of  the 
nonlinear  material  model  in  the  continuum  mechanics  based  ANCF  shell  element  developed  in 
this  investigation.  The  Lame  constants  p  and  X  are  assumed  to  be  8.077x10^  Pa  and 

1.212x10^  Pa,  respectively.  The  corner  of  the  quarter  cylinder  is  connected  to  the  ground  by  a 
spherical  joint.  The  deformed  shapes  of  the  quarter  cylinder  under  the  effect  of  gravity  are  shown 
in  Fig.  1 1 ,  in  which  the  shell  experiences  large  deformation.  The  global  X,  Y  and  Z  position  at 
point  A  shown  in  Fig.  11  are  presented  in  Figs.  12-14  for  different  number  of  elements,  and  it  is 
shown  that  use  of  32x32  elements  leads  to  the  convergent  solution. 

6.  SUMMARY  AND  CONCLUSIONS 

In  this  investigation,  a  continuum  mechanics  based  bi-linear  shear  deformable  shell  element  is 
developed  using  the  absolute  nodal  coordinate  formulation  for  the  large  deformation  analysis  of 
multibody  shell  structures.  The  elastic  forces  are  formulated  using  the  continuum  mechanics 
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approach  which  allows  for  the  consideration  of  nonlinear  material  models  such  as  hyperelasticity 
in  the  shell  element  in  a  straight  forward  manner.  By  doing  so,  the  stress  distribution  along  the 
thickness  is  considered  and  the  elastic  forces  are  evaluated  as  a  continuum  volume.  It  is 
demonstrated  by  some  numerical  examples  that  use  of  the  combined  EAS  and  ANS  approach 
alleviates  the  thickness  locking  effectively  exhibited  in  the  continuum  mechanics  based  bi-linear 
shear  deformable  ANCF  shell  element,  while  the  transverse  shear  locking  as  well  as  the  in-plane 
shear/normal  strain  locking  can  be  alleviated  by  the  ANS  and  EAS,  respectively.  The  locking- 
free  shear  deformable  element  developed  in  this  investigation  can  be  used  for  modeling  large 
deformable  multibody  shell  structures  in  challenging  engineering  problems  with  material 
nonlinearities  that  include  the  modeling  of  rubber  tires  for  vehicle/terrain  interaction,  wind 
turbine  rotor  blades  made  of  composite  materials,  and  biomechanics  applications  with 
biomaterials. 
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TABLE  CAPTIONS 

Table  1.  Models  for  continuum  mechanics  assumption 

Table  2.  First  ten  dimensionless  eigenfrequencies  £1  =  co/a>Qof  the  free  (FFFF)  square  plate 

FIGURE  CAPTIONS 

Figure  1 .  Kinematics  of  bi-linear  ANCF  element 
Figure  2.  Sampling  points  for  assumed  natural  strain 

Figure  3.  Deformed  shape  of  a  cantilevered  plate  subjected  to  a  large  transverse  tip  load 
Figure  4.  Numerical  convergence  with  large  deformation  (initially  flat) 

Figure  5.  Deformed  shape  of  a  cantilevered  plate  subjected  to  a  large  transverse  tip  load 
Figure  6.  Numerical  convergence  with  large  deformation  (initially  curved) 

Figure  7.  Deformed  shape  of  pinched  semi -cylindrical  shell 

Figure  8.  Load-deflection  curve  of  pinched  semi-cylindrical  shell 

Figure  9.  Deformed  shape  of  slit  annular  plate  subjected  to  lifting  force 

Figure  10.  Load-deflection  curve  of  slit  annular  plate  subjected  to  lifting  force 

Figure  1 1 .  Deformed  shapes  of  quarter-cylindrical  shell  pendulum 

Figure  12.  Global  X-position  at  Point  A 

Figure  13.  Global  Y-position  at  Point  A 

Figure  14.  Global  Z-position  at  Point  A 
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Table  1.  Strain  modification  for  the  continuum  mechanics  based  ANCF  bi-linear  shell  element 


Model 

name 

ANS  for  transverse 
shear  strains 

.ANS  ^ANSx 

HAS  for  in-plane 
strains 

/^eas  has  ^EASx 

ANS  for 

transverse 

normal  strain 

(^r) 

HAS  for 

transverse 

normal  strain 

(C) 

Combined 
EAS  and  EAS 
for  transverse 
normal  strain 

/  ANS/EAS  X 

1  zz  ^ 

Model- 1 

- 

- 

- 

- 

- 

Model-2 

Y 

- 

- 

- 

- 

Model-3 

Y 

- 

Y 

- 

- 

Model-4 

Y 

Y 

Y 

- 

- 

Model-5 

Y 

Y 

- 

Y 

- 

Model-6 

Y 

Y 

- 

- 

Y 
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Table  2.  First  ten  dimensionless  eigenfrequencies  O  =  (o/Oq  of  the  free  (FFFF)  square  plate 


1x1 

2x2 

4x4 

8x8 

16x16 

Analytical 

solution 

1.4328 

1.4383 

1.4014 

1.3734 

1.3548 

1.3646 

101.5117 

2.3477 

2.1142 

2.0141 

1.9721 

1.9855 

101.7252 

3.1993 

2.7186 

2.5169 

2.4415 

2.4591 

101.7288 

4.0870 

3.8097 

3.6023 

3.5267 

3.5261 

115.9849 

4.0906 

3.8172 

3.6153 

3.5524 

3.5261 

115.9849 

7.5903 

7.4061 

6.7273 

6.2849 

6.1900 

138.3927 

98.8943 

8.4204 

6.7365 

6.2862 

6.1900 

138.6304 

98.8949 

8.4332 

6.7394 

6.5196 

6.4528 

152.2271 

101.7244 

8.7211 

7.4840 

7.1160 

7.0181 

152.4032 

101.7288 

10.5372 

8.5216 

7.9556 

7.8191 
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Figure  1.  Kinematics  of  bi-linear  shear  deformable  ANCF  shell  element 
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Figure  2.  Sampling  points  for  assumed  natural  strain 
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Figure  3.  Deformed  shape  of  a  cantilevered  plate  subjected  to  a  tip  force 
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Figure  4.  Numerical  convergence  of  ANCF  shell  solutions:  cantilevered  plate  subjected  to  a  tip 

force  with  large  deformation 
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Figure  5.  Deformed  shape  of  a  quarter  cylinder  subjected  to  a  tip  force 
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Figure  6.  Numerical  convergence  of  ANCF  shell  solutions:  quarter  cylinder  subjected  to  a  tip 

force 
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Z(m) 


Figure  7.  Deformed  shape  of  pinched  semi-cylindrical  shell 
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Deflection  (m) 


Figure  8.  Load-deflection  curve  of  pinched  semi-cylindrical  shell 
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Z(m) 


Figure  9.  Deformed  shape  of  slit  annular  plate  subjected  to  lifting  force 
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Deflection  (m) 


Figure  10.  Load-deflection  curve  of  slit  annular  plate  subjected  to  lifting  force 
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Point  A 


Figure  11.  Deformed  shapes  of  quarter-cylindrical  shell  pendulum 
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Figure  12.  Global  X-position  at  Point  A 


38 


Global  Y-position  (m) 


UNCLASSIFIED:  Distribution  Statement  A.  Approved  for  public  release.  #24515 


Figure  13.  Global  Y-position  at  Point  A 
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Figure  14.  Global  Z-position  at  Point  A 
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